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Abstract 


This paper presents an overview of the theories and computer implementation aspects 
of phase field models (PFM) of fracture. The advantage of PFM over discontinuous ap- 
proaches to fracture is that PFM can elegantly simulate complicated fracture processes in- 
cluding fracture initiation, propagation, coalescence, and branching by using only a scalar 
field, the phase field. In addition, fracture is a natural outcome of the simulation and 
obtained through the solution of an additional differential equation related to the phase 
field. No extra fracture criteria are needed and an explicit representation of a crack surface 
as well as complex track crack procedures are avoided in PFM for fracture, which in turn 
dramatically facilitates the implementation. The PFM is thermodynamically consistent 
and can be easily extended to multi-physics problem by ’changing’ the energy functional 
accordingly. Besides an overview of different PFMs, we also present comparative numerical 


benchmark examples to show the capability of PFMs. 
Keywords: Phase field, Brittle fracture, Computer implementation, Finite element method, 
Hydraulic fracture 
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1 Introduction 


The prediction of material failure is of major importance in engineering and material science 
(Rabezuk 2013). Consequently, many numerical approaches have been developed to handle 
fracture problems in recent years. They can be classified into two categories namely discon- 
tinuous and continuous approaches. Discontinuous approaches introduce a strong discontinuity 
in the displacement field. Typical discontinuous methods are discrete crack models (Ingraffea 
and Saouma 1985), the extended finite element method (XFEM) (Moés and Belytschko 2002; 
Moés et al. 1999), generalized finite element method (GFEM) (Fries and Belytschko 2010), and 
the phantom-node method (Chau-Dinh et al. 2012; Rabczuk et al. 2008). Most discontinuous 
approaches to fracture require a method to represent the crack’s topology — such as triangular 
facets (Zhuang et al. 2012) or level sets (Zhuang et al. 2011) — and associated crack tracking 
algorithms. Obtaining complex fracture patterns such as crack branching and crack interactions 
require additional criteria (the TLS approach (Moés et al. 2011) is a discontinuous approach that 
does not require special criteria for branching). Some of the above methods also employ special 
schemes to ‘treat’ crack tip singularity, which improves accuracy but imposes other difficulties 
such as numerical integration. 

In contrary to the discontinuous approaches, the continuous approaches to fracture do not 
introduce discontinuities in the displacement field. Popular continuous approaches include gra- 
dient damage models (Peerlings et al. 1996), screened-Poisson models (Areias et al. 2016a,b) 
and phase field models (PFMs) (Borden et al. 2012; Miehe et al. 2010a,b). All these models in- 
troduce an intrinsic length scale and smear the fracture over a localization band of finite width. 
This paper focuses on the phase field model of brittle fracture. Phase-field models of brittle 
fracture can be traced back to the late 1990s and received extensive development in theory and 
computer implementation. The crack is diffusively represented by a scalar field (phase field) and 
the evolution equation of the phase field is used to model crack propagation. 

The phase field models recently received extensive attention because they can elegantly 
simulate complicated fracture processes including crack initiation, propagation, coalescence, and 
branching quite naturally. In addition, the fracture evolution can be simulated on a fixed mesh. 
The phase field model avoids the laborious task to track crack surface, which is especially tedious 
in 3D. 

In the past decade, researchers made a huge effort to develop novel, efficient, and accurate 
phase field models and achieved enormous progress. This paper aims to review the advances in 
this direction. It is organized as follows. Section 2 gives an overview of the theories used in the 
phase field modelling. Section 3 is devoted to the different discretization methods in fracture 
modelling. Section 4 gives detailed FE discretization and some details about the computer 
implementation in phase field modeling. An overview of the state-of-the-art applications and 


extensions of the phase field models is given in Section 5 followed by Section 6 which presents 


some representative numerical examples to show the capability and practicability of PFM in 2D 


and 3D. Finally, we end with concluding remarks and future directions in Section 7. 


2 Theories of phase field models for fracture 


In this section, we provide the basic theories of the phase field method for fractures. The 
phase field models are different in the ‘physics’ and ‘mechanics’ communities. The difference 
between these two models lies in that in the physics community, the models commonly come 
from the Landau-Ginzburg phase transition (Aranson et al. 2000) without introducing the idea 
of length scale to diffuse the crack. In addition, the free energy in the models don’t contain the 
fracture energy. The PFMs in the physics community are often used to model dynamic fractures 
(Aranson et al. 2000; Henry and Levine 2004; Karma et al. 2001). However, the phase field 
models in the mechanics community (Miehe et al. 2010a,b), can be regarded as the extension of 
Griffith's fracture theory although the crack evolution is similar to the Landau-Ginzburg phase 
transition while a clear length scale parameter is used. The free energy used for variational 


analysis contains the fracture energy, which is regularized by using the length scale parameter. 


2.1 Physical models based on Landau-Ginzburg phase transition 
2.1.1 Aranson, Kalastky, Vinokur model, 2000 


The fracture model presented by Aranson et al. (2000) is among the first phase-field like descrip- 
tions of crack propagation in brittle materials. Their model focuses on Mode-I fracture and they 
implemented their approach in 2D. The displacement u satisfies the standard elastodynamic 
equation with a damping term: 

pü = nù +V -0 (1) 


where p is the density of material, V and A are the divergence and Laplace operators, ø is the 
stress tensor, y > 0 is a viscous damping parameter, ù = ĝðu/ðt, and ü = 0%u/0t* with t being 
the time. 

Aranson et al. (2000) defined a local-order parameter s (a field parameter): s = 1 outside 
the crack (no defects) and s = 0 inside the crack (all the atomic bonds are broken). The order 
parameter s is assumed much larger than the inter-atomic distance, thereby obeying the contin- 
uum description of fracture. Materials cannot bear tensile stresses and fail at an order parameter 
below the critical value se (Aranson et al. 2000). The stress-strain relation subsequently differs 
from that used in the isotropic linear elasticity o = C : e by introducing the dependency on the 
order parameter s: 

o=sC:e+ xslI (2) 


where C is fourth-order elasticity tensor of material, e is the strain tensor, I is the second- 
order identity tensor, and x > 0 is an additional material parameter describing the hydrostatic 
pressure due to creation of new defects. 

The order parameter s is assumed to be governed by pure dissipative dynamics and $ = 
—dE/ds with E being a "free-energy” type functional. Based on the Landau's phase transitions 
(Landau et al. 1980), the simplest form of E reads 


B= [ [P(6) DVsf| do (3) 


where P is a polynomial functional and D, > 0 is an adaptive constant. 
In the fracture model of Aranson et al. (2000), the order parameter evolution equation is 


naturally expressed as 
§ = D,As — s(1— s) [a (1 + (tr(e) — b)s) — cu - Vs] (4) 


where a, b, and c are model parameters. 

In general, the Aranson, Kalastky, Vinokur model can simulate multiple fracture behaviors 
such as crack initiation, propagation, branching, dynamic fracture instability, sound emission 
and fragmentation. However, some discrepancies still exist between the numerical predications 
and experimental observations of Mode I crack propagation in a rectangular strip of finite width 
(Aranson et al. 2000). 


2.1.2 Karma, Kessler, Levine model, 2001 


Karma et al. (2001) proposed a phase-field model for Mode III (antiplane shear) fracture. In their 
model, only the out-of-plane displacement component exists in a 2D setting. Their approach 
for fracture emanates from the original phase field models for solidification and the basic free 


energy function is expressed as (Hakim and Karma 2009; Karma et al. 2001) 
1 
Ets) = f lat (Vole) — V.)  V(s) + 2D,Vsl" | 40 (5) 
Q 


where g(s) is a function of the order parameter s and satisfies g(s) > 0 for 0 < s < 1. V(s) = 
s?(1 — s?)/4 is the so-called Ginzburg-Landau double-well potential (Ambati et al. 2015b). In 
addition, V, denotes the strain energy threshold for crack initiation, and D, > 0 is the positive 
constant identical to that in Eq. (3). Note that the function g(s) = s?** is used in Karma 
et al. (2001) with a4 a dimensionless coefficient. Subsequently, the energy functional is used 
to governing the main equations of systems, including the momentum balance equation, the 


stress-strain relation and the phase-field evolution equation. 


That is, the variation of the functional E with respect to u results in the equilibrium equation: 


pu = V - 0. The stress-strain form is subsequently represented as follows: 


c (u, 5) = g(s) - g(s)C: e (6) 
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'The proposed fracture evolution law reads 75 — ae with 7 > 0 being a kinetic modulus. 


The nonlinear evolution equation is thus expressed as follows: 
Ts = D,As — V'(s) — g'(s) (Wo(e) — V.) (T) 


2.1.3 Henry and Levine model, 2004 


Henry and Levine (2004) subsequently extended the original model of Karma et al. (2001). A 
modification of the elastic energy is seen in Henry and Levine (2004) for better simulating the 
crack growth of Mode I and II under 2D plane strain condition. The total energy functional 
retains the basic structure of the Karma, Kessler, Levine model (KKL model) and still prevents 


the compressed region of the elastic body from cracking. The total energy is expressed as follows, 
- 1 
E(u,s) — J (9) (He) — v.) +V(s)+ 5 D.lVsl dQ (8) 
Q 


where W(e) is the modified elastic strain energy and it is identical to the standard elastic energy 
density for a positive volume strain. However, for a negative volume strain, a breaking symmetry 


term is introduced: 


Yole) if tr(e)>0 
e= Y aste) ¿0 Kt(e) ie) M 


where Wo(e) = 3Ae7; + we;,, A and u are Lamé constants, K = (A + u)/2, and az > 1 is an 


arbitrary coefficient. It should be noted that the coupling function g(s) used in the model of 
Henry and Levine (2004) is different from that in Karma et al. (2001) and here g(s) = (4—3s)s?. 


The stress-strain relation used in Henry and Levine model is 


i 1 
o(u, 5) = 9(s) (10) 
while the evolution equation of s is modified as 
T$ = D,As — V'(s) — g'(s) (vo) = v.) (11) 


In summary, the Henry and Levine model can accurately reproduce different behaviors of 
cracks such as branching, experimentally observed oscillating cracks, and well observed super- 
critical Hopf bifurcation (Henry and Levine 2004). 


2.2 Mechanical models based on Griffith’s fracture theory 
2.2.1 Variational approach to fracture 


The original phase field models (Miehe et al. 2010a) in the mechanics community is developed in- 
dependently for quasi-static fractures. These approaches originate from the variational approach 
of brittle fracture proposed by Francfort and Marigo (1998) and refer to the regularization for- 
mulation used by Bourdin et al. (2000). 

Similar to the models in the physics community, the variational approach to fracture also 
requires constructing an energy functional that governs the entire fracture process. That is, crack 
initiation, propagation, and branching is a minimization of the energy functional (Francfort and 
Marigo 1998): 


bet) = / Wole)dQ +G. / ds (12) 


where G. is the fracture toughness and also referred to as the critical energy release rate. IT € Q 
is an admissible crack set and the displacement field u is discontinuous across I’. 

To make the variational formulation amenable to a numerical implementation, Bourdin et al. 
(2000) proposed a regularized version of the variational formulation. Without treating the free 
discontinuity sets of the displacement field u, an auxiliary scalar s(a,t) (a the position vector) 
is employed to diffusely represent the sharp fracture geometry. Identical to the physical phase 
field models, the scalar field s continuously transits between s = 1 (intact material) and s = 0 
(fully damaged material). The main advantage of introduction of the phase-field is that the 
representation of cracks is no longer mesh or geometry based (Kuhn 2013). The regularized 
approach can be easily implemented in a finite element framework and retain the main advantage 
of the variational formulation to model cracks. In the regularization formulation of Bourdin et al. 
(2000), the modified energy functional reads 


1 
E(u, s) = [eo + n)WVo(e(w))dQ + c. | P — s)’ + eV s|? dQ (13) 

Q Q LHE 
where the parameter 0 < yn < 1 is used to avoid numerical singularity when the material is 
broken. Another parameter e has the dimension of a length and controls the transition zone 
between the fully damaged and intact body. To obtain the displacement u and phase field s, 
minimization of the modified energy functional is also required. It should be noted that the 


regularized approach can be recovered to the original variational approach to fracture when 


e — 0 in the sense of I'-convergence (Braides 2006). 


2.2.2 Kuhn and Müller model, 2008 


Kuhn and Müller (2008) reinterpreted the crack variable as a phase field order parameter and 
they also regarded cracking as a phase transition problem. By applying the thermodynamics 
framework of order parameter based models, Kuhn and Müller (2008) reformulated the mini- 


mization problem in Eq. (13) and used the stress equilibrium equation V -ø = 0 with 


Yole) 
ðe 


o — (s? +n) = (+ n)C: e (14) 


The crack propagation is simulated by using the evolution equation of the order parameter 


s. The Ginzburg-Landau type evolution equation (Kuhn and Müller 2008) then reads 


pec ESO ES (zeas pc x 3l (15) 


Equation (15) is used under the irreversibility constraint of crack evolution and with a mo- 
bility parameter M. M > 0 governs the dissipation of stable crack propagation. At a finite 
value of M, the crack model can be also regarded as a viscous quasi-static model (Miehe et al. 


2010a) and in a limit case M — oo quasi-static crack propagation thereby obeys 


] — 
2sWo(e) — Ge (zeas += : =0 (16) 


In Kuhn and Miiller (2008), the numerical model is implemented within the finite element 
framework and an implicit Euler scheme is used for time integration. Note that the Kuhn and 
Müller model (Kuhn and Müller 2008) has the similar forms and derivations of the Karma, 
Kessler, Levine model (Karma et al. 2001) and Henry and Levine model (Henry and Levine 
2004). 


2.2.3 Amor, Marigo, Maurini model, 2009 


The originally developed phase field model (e.g. Eq. (13)) does not distinguish between fractures 
due to compression and tension. Therefore, some unrealistic crack patterns are reported in 
Bourdin et al. (2000). To avoid unrealistic simulations, elastic energy must be decomposed 
into different parts. Amor et al. (2009) made the preliminary contribution to prevent phase 
field representation of cracks due to compression. Amor et al. (2009) modified the regularized 


formulation of Eq. (13) by introducing an additional decomposition of the elastic energy Wo(e) 


into volumetric and deviatoric parts. That is, Yo = Vj + Y3 with 


WE (e) = 5 Kaltr(e))3 + ule? : e?) an 


WE (e) = 5 (tr()) 


where Ka = A+ u, €? = e — ztr(e)I, and operators (*), = 5(* |» |). Amor et al. (2009) 


subsequently proposed the modified energy functional as 


1 
E(u,s) = f ((s? +n) Ug (e) + Wp (&)) dQ + Ge | P — sy? +¢|Vsl|?| dO (18) 
Q Q 
Amor et al. (2009) applied an alternate minimization algorithm to solve for the energy 
functional (18) and achieved local minimization. They solved a series of minimization sub- 
problems on u at a fixed s, and vice versa on s at a fixed u until convergence. Ambati et al. 
(2015b) pointed out that a good outcome of the energy split of Vo is the resulting stress-strain 


relation: 
OW; " OW, 

Oe Oe (19) 
=(s? +1) [Knltr(e)) E + 2ue?] + Kn(tr(e))- 


o(u, s) =(s? + n) 


However, the phase evolution equation suggested by Ambati et al. (2015b) is not the vari- 
ational outcome of the energy functional (18). Instead, the Ginzburg-Landau type equation is 
still used and Wo in Eq. (15) is replaced by Vj: 


¿=-M ESO -Q, (zeas qe x 3l (20) 


Thus, the evolution of s in Eq. (20) is driven only by the positive parts of the elastic energy. 
That is, only volumetric dilatation and shear deformation can produce cracks. It should be 
noted here that the Amor, Marigo, Maurini model does not obey the P'-convergence although 
the energy decomposition model (Amor et al. 2009) produces adequate numerical results. The 
reason is that it is not clear what kind of functional and physical process are to be recovered 


when € > oo. 


2.2.4 Miehe et al. model, 2010 


Miehe et al. (2010a,b) presented a thermodynamically consistent phase field model of brittle 
fracture. The model of Miehe et al. (2010a,b) is another important development after the 
variational framework of fracture was presented by Francfort and Marigo (1998) and Bourdin 
et al. (2000). This model is based on continuum mechanics and uses an auxiliary scalar field 
@ € [0,1] (the so-called phase field). The additional phase field is used to smear the sharp 


10 


crack shape. We can also regard $ = 1 — s here with relation to the former phase field models 
such as the Amor, Marigo, Maurini model. ¢ = 0 and Y = 1 represent the intact and fully 
broken states, respectively. A length scale parameter ly is used in Miehe et al. model and ly 
controls the transition region between the fully broken and intact bodies. Another elastic energy 
decomposition method is used with Vo = Vj + W7 based on the spectral decomposition of strain 
tensor e: e = 35 ,emi & ny. {e}}_, and {n}3_ are the principal strains and their directions. 


'The decomposed elastic energy reads 


WE (e) = 5A (tr) ute (ed) (21) 


where €, = (enn; @ ny. 


The modified elastic energy V, used in the total energy functional is expressed as follows, 


ve = [0 — 4)? m) Vj + V (22) 


The governing equation of the displacement is V -øo = 0 and the stress tensor is expressed 


as 
OV; | OW, 
- 2 0 0 
The evolution equation of phase field in the model of Miehe et al. (2010a,b) is 
P Ge , 
Ep =2(1— p)W (E) + — (6 — Ag) (24) 


lo 


where £ is a viscosity parameter. 

It should be noted that if lọ = 2e and € = 3, the structure of Eqs. (23) and (24) resembles 
those equations in the Amor, Marigo, Maurini model (Amor et al. 2009). The energy decomposi- 
tion of Miehe et al. (2010a,b) is different from that of Amor et al. (2009). Miehe et al. (2010a,b) 
only distinguishes compressive and tensile parts and imposed a degradation of the compressive 
energy component. However, the deviatoric stress is considered in Amor et al. (2009). Another 
feature of Miehe’s model is that Eq. (23) results in a highly non-linear stress-strain relation- 
ship. Therefore, a much higher computational cost is needed compared with former physical 
and mechanical phase field methods. In order to ensure a monotonically increasing phase-field, 
the irreversibility condition must be imposed as a constraint during compression or unloading. 
Miehe et al. (2010a,b) proposed a new approach where a history-field H is defined in a loading 
process: 

H (z,t) = max yf (e(x,7)) (25) 


TE(0,t] 
The introduction of the history field enhances the phase-field formulations and overcomes 


some implementation difficulties. Replacing Yọ in Eq. (24) by H, the finally used evolution 
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equation reads E 
E$ = 2(1— 9)H + p» (ó — I5Ad) (26) 


The common isotropic, anisotropic, and hybrid phase field models in the mechanics commu- 
nity are summarized in Table 1. Thus, the differences in governing equations and driving forces 


for different PFMs can be easily identified by readers. 


2.3 Phase field approximation of the sharp crack topology 
2.3.1 Classic crack surface density function 


For the crack phase field s(x,t) or ¢(a,t), let us follow Miehe et al. (2010a,b) and introduce a 
regularized functional A(@) or A(s) such that 


Aal) = f yd, VOA ~ l TUE / dS — A, 


(27) 
Aa(s) = [s van e f san - ] 35-4. 
Q Q r 
where 7(¢, Vo) or y(s, Vs) is the crack surface density function, and A, is the crack area. 
The crack surface density function approximates the Dirac-delta 6, along the crack I, and is 
composed of the crack phase-field and its spatial gradient. 
Bourdin et al. (2000) proposed the original form of the crack surface density function as 


y(s, Vs) = i (1— s)? + e| Vs? (28) 


de 
Miehe et al. (2010a,b) followed the original form of Bourdin et al. (2000) but used different 


parameters and phase field: 
1 l 
1(¢,V9) = = # + SIVOP? (29) 
2lo 2 


In the physical phase field models, the crack surface density function follows a Ginzburg- 
Landau double-well potential (Hakim and Karma 2009; Karma et al. 2001): 


s?(1— s?) 
4 


1 
qls, Vs) = + 5PslVsl’ (30) 


2.3.2 Generic geometric crack function 


Kuhn, Schliiter, Miiller model, 2015 Kuhn et al. (2015) proposed a generic energy density 


functional V of a phase field fracture model: 


V = (g(s) + 1) Vo 4 E e | dv) (31) 
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Table 1: Different forms of phase field methods in mechanics community 


Type 


Isotropic 


Anisotropic (Miche 
et al. 2010a) 


Anisotropic (Amor 
et al. 2009) 


Hybrid (Ambati et al. 


2015b) 


Hybrid (Zhang et al. 
2017) 


Hybrid (Zhou et al. 
2019a) 


Governing equations 


V-o0+b=0 
2lo 
Gc 


BAd+o= —(1-¢)H 


o(u,¢) = (1 - p)? 
Oe 


_ oe 
^ Oe 
V.o-b-0 
20H 202,  2loH 
+ 1 lV 
(E +1) e- ive = 70: 
Out Ove 
= (1 2 € E 
o(u,6)=1- y SE 4 
dat do 
D= 
ðe ðe 
V-o+b=0 
2l 
BA +o = c0 -9)H 
Ovg , Oba 
=(1 29V9 , 0 
olup) = (1-4)? 2 + 8 
LED. OE 
^. 0€ ðe 
V-o+b=0 
2 2lo 
bAptp= 4 $)H 


o(u,¢) = (0 - oe 
Oe 


_ Oc 
~ 8e 
V.o+b= 
Hi, Ha ho cds 
a ES cs.) ag "id i 
a 
0(u, 0) = 1 - 99 8 
p= 9€ 
Oe 
V:-o0+b=0 
2l 
BAd+o= c0 - 9)H 
a 
o(w,4) = (1-9)? 2 
€ 
_ a0 
^ Oe 
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Definition 


H = max Yo (e(#,7)) with 
TE[O,t] 


bo = 3 (tr(e))? + utr(e?) 


H = max V (e(#,7)), 
T € [0,1] 


+— 2005 ¿-_ Ove 
ce (1 $) ðe y E ðe 3 


be (e) = }(tr(e))} + utr (e3) 


dug __ æ 
ot =(1 py? 52-0 = ac 


Vg = 3 Kn ((tr(e))4 + ule : ede), 
Vg = 4 Kn ((tr(e))? with Kn a 
physical parameter 


H = max wt (e(a,7)), 
T € [0,1] 


Yo = A(tr(e))? + ptr(e2), 
vi (e) = A(tr(e))2 + ptr (e2) 


Hi = E A(tr(&))2., 


H = t 2], Ger and Ge 
2 mex nerie] 1 and Gerr 


are the critical energy release rate of 
modes I and II 


H= max v» (e(@,7)), bp = 


3 2 1 /H(5ip—£jp) 
Pj y» 2G ( cose H [A(E1p4 


HE2p + €3p) + M(Eip + Ejp)] tang — c)? 


where g(s) is the energetic degradation function with g(1) = 1 and g(0) = 0. The function w(s) 
models the local fracture energy and cw is a normalization constant that results in the integral of 


(42 + €|Vs|?) /(2c,,) over the fractured domain converges to the surface measure of the crack 
set when e — 0 (Kuhn et al. 2015). Kuhn et al. (2015) concluded two types of w(s), namely, the 
double well function w(s) = 16s?(1 — s?) and monotonous functions w(s) = (14- 8s)(1— s) with 
B € [-1,1]. The convex quadratic function with 6 = —1 is mostly used, such as in the model of 
Miehe et al. (2010a). In addition, models with @ = 0 can be found in Hossain et al. (2014). 
Kuhn (2013) compared the function of w(s) = 16s?(1— s?) and w(s) = (1— s)?. The function 
w(s) of these two forms is shown in Fig. 1. Their observations proved that the double well 
potential naturally models the irreversibility of fracture processes to a certain extend because 
of local maximum acting as an energy barrier between the broken and undamaged phase. In 
contrast, the crack phase field needs additional irreversibility to prevent crack healing for w(s) — 


(1 — s)? due to lack of energy barrier between the broken and undamaged states. 


—— w(s) 2 s 
—---w(s) = 16s (1-5) | 


Fig. 1 Different degradation functions 


Wu's model, 2017 Wu (2017) proposed another generic form for the crack surface density 
function (o, Vo): 


1(0, vo) = = (ZrO uve?) (32) 
with 


ed / Vrae (33) 


where &(ó) is the geometric function that characterizes homogeneous evolution of the crack 
phase-field, b is an internal length scale parameter that regularizes the sharp crack, cy > 0 is 


a scaling parameter by which the regularized functional Ag(@) can be recovered to the crack 
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surface A, when b — 0. The proposed crack geometric function «(¢@) must obey 


K(0)=0, and k(1)=1 (34) 


2.4 Energetic degradation function 
2.4.1 Classic forms 


As aforementioned, g(s) is a degradation function that models the degradation of the elastic 
energy when a fracture initiates. Seen from a simple stress-strain relation ø = (g(s) + n)e, the 
function g(s) also models the decrease in the stiffness in a broken material when the phase field 
s — 0. In addition, the energetic degradation function g(s) must monotonically increase with 
g(1) = 1 and g(0) = 0. Because the derivative g'(s) enters the evolution equation of the phase 
field, g'(0) = 0 must be imposed to eliminate the elastic driving term ae = Lg (s)e : e when 
s=0. 

Some of the classic degradation functions have been mentioned in the overview of the most 
important phase field models such as the quadratic function g(s) = s? in Bourdin et al. (2000), 
and the quartic function g(s) = 4s? — 3s* in the KKL model (Karma et al. 2001). In addition, 
Kuhn et al. (2015) adopted a cubic function: 


g(s) = 3s? — 25? (35) 


2.4.2 Wu's model, 2017 


Wu (2017) proposed a unified phase field damage model for brittle fracture. The elastic energy 
density V is modeled as 
V = v(ó)Vo (36) 


The monotonically decreasing energetic function v(p) describes degradation of the initial 
energy Vo with the crack phase-field evolution. The function v($) obeys (Miehe et al. 2010a) 


v (6) « 0, and #(0)=1, «v(1)=0, v(1)=0 (37) 


The generic form of the energetic degradation function proposed by Wu (2017) is given by 


(1 9)" 
T= 6 QU) 


v(p) = (38) 


with the exponent pı > 0 and continuous function Q(o) > 0. 


15 


For a strictly positive function Q(¢), the following polynomial is considered 


Q(¢) = a1¢ + aiaa? + a102030? +--+ = o40P(Q) (39) 
P(9) = 1+ a2¢ + 09030" +--+ (40) 
More details about the determination of the coefficients o4, o», --- can be referred to Wu 


(2017). 


2.4.3 Sargado et al. model, 2017 


Sargado et al. (2017) proposed a new parametric family of degradation functions to increase the 
accuracy of phase field models in predicting critical loads. Expect a better prediction of the 
crack initiation and propagation, their additional goal is to preserve linear elastic response in 
the bulk material prior to fracture. Their numerical examples indicated the superiority of the 
proposed family of functions to the classical quadratic degradation function. 

The family of degradation functions proposed by Sargado et al. (2017) are based on three 


parameters: 
1 — ete)" 
gel) = (1 — w) —— + wf.(9) (41) 


where k, n and w are parameters and k > 0, n > 2 and w € [0,1]. The function fe acts as a 


corrector term. If w = 0, the proposed function has only two parameters and follows (a) g.(@) 


is monotonically decreasing, (b) ge(0) = 1, ge(1)=0, and (c) g} < 0, gL(1) = 0. 


2.5 On constitutive assumptions 


As aforementioned, the model of Amor et al. (2009) and the model of Miehe et al. (2010a,b) 
achieve variationally consistent constitutive model for stress. In addition, the crack driving force 


D, in the model of Amor et al. (2009) is derived by consistent variation: 
1 2 Bis, up 
D, = 5 Ku(tr(e))3 + ule? : e?) (42) 
Consequently, only positive volume changes and distortion of the shape can produce crack 
evolution. For the model of Miehe et al. (2010a,b), the crack driving force D;: 


D, = 5Mtn(e))2. uta(ed) (43) 


Both volumetric-deviatoric splits (Amor et al. 2009) and spectral decomposition (Miehe 
et al. 2010a,b) fully preserve the variational character of the phase field method. The en- 


ergy split scheme ensures the energy contribution in the evolution of phase field and the non- 
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interpenetration condition on crack surfaces. This means that the normal stress on the crack 
should be non-positive and the shear stresses along a frictionless crack should vanish (Strobl and 
Seelig 2016). On this, some variationally non-consistent constitutive equations for the stress are 
developed. A constitutive model based on a crack orientation dependent decomposition of ø is 
mentioned in Strobl and Seelig (2016) where e, and op represent the stress affected and unaf- 


fected by the phase field respectively. The stress decomposition depends on the normal strain 


Enn: 


fo = g(s)(Atr(&)1 + 2ue), if Enn > 0 (44) 


c, = g(s)(Atr(e)1 + 2ue) + (1 — g(s))(A--2u)(e : NIN, if Em > 0 
where N =n, 9 n, with n, being the normal vector of the crack surface. 
A more sophisticated stress model can be seen in Strobl and Seelig (2015) and only the 
stiffness normal to the crack is degraded. Tensile and shear stresses vanish on the fully developed 


crack surfaces. This results in the following stress: 


2 


A+ 2u 


2 


) tr(e)1 + 2ue + (g(s) — 1) E t > À ) (tr(e)N + (e : N)1) 


on = (à+ (0) = x 


2 


A+ 2u 


+ 4(1 — g(s)) (2+2- ) (e: NIN + u(g(s) 1N -e+e: IN) (45) 


For compression, only shear stresses parallel to the crack surfaces are degraded. The modified 


tensor is expressed as 
c, = Atr(e)1 + 2ue + 4u(1 — g(s))(e: NN + u(g(s) - 1) (N -e+e- N) (46) 


Equation (46) shows a transverse isotropy with symmetry about the crack normal and the 


stiffness degradation is described by A, u and the phase field. 


2.6 Governing equations of PFM in strong form 


In summary, the quasi-static phase-field models, which are developed from the regularized vari- 
ational formulation of Bourdin et al. (2000), can be classified into three basic types, namely, the 
isotropic, anisotropic, and hybrid (isotropic-anisotropic) types. The isotropic and anisotropic 
formulations are the natural results of the variational principle and the basic structures can be 
expressed as follows. 


(1) Isotropic formulation 


V.oc-b-0 


—¢Ab+¢= (1 —d)H 


c 
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where o (u, ġ) = (1 — oem and H = max Yo (€(#,T)). 


TE(0,t] 
(2) Anisotropic formulation 


V-o+b=0 


48 
- Roao+o= 20 =d)H" x 


where a(u,ó) = (1 — 0) E + GE and H* = max dij (e(as 7). 
TEO 


It should be noted again that the isotropic formulation is easier to implement and has lower 


computational cost than the anisotropic formulation because the stress-strain constitutive re- 
lationship is linear in the isotropic approach. However, non-linear stress-strain relation exists 
in the anisotropic formulation. Despite of the lower implementation cost, the isotropic formu- 
lation can be used only in some very simple problems because fractures in compression and 
interpenetration of the crack faces are allowed. The anisotropic formulation which uses energy 
split technology, naturally overcomes the drawbacks emanating from the isotropic formulation. 
However, because of non-linearity from strain decomposition, the numerical implementation of 
the anisotropic models is laborious. 

Ambati et al. (2015b) developed the hybrid phase field method, which retains a linear mo- 
mentum balance equation. Thus, a favorable computational cost of the isotropic methods is 
retained. However, the evolution equation of the anisotropic methods is also retained to avoid 


fractures in compression. The hybrid model of Ambati et al. (2015b) reads 


V.a+b=0 
— eA + ¢ = (1 dm 
m (49) 


= (1-41 
olup) - (1 - o 
Va: V cWy—60-0 
Wu (2017) also used the hybrid formulation but the driving force of the crack field is modified 


as 
1 


D, — 2E0 Da 


where Ey is the undamaged Young's modulus and Geg is the equivalent stress. More detailed 


(50) 


descriptions about the equivalent stress can be seen in Wu (2017). 

Another variationally non-consistent phase field method can be seen in Zhang et al. (2017). 
The driving term A/G. in the evolution equation is decomposed into two parts: H,/G.; and 
Hə/Gerr with 

H, = max A(tr(e))% (51) 


TE(0,t] 
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and 
H, = max jutr[(e)?] (52) 


TE[0,t] 
Zhang et al. (2017) regarded the parameters Ger and Gerr as the respective strain energy 
release rates due to mode I and mode II deformations at a crack tip. However, their definition 
cannot practically distinguish between the tension and shear fractures. The governing equations 


of Zhang et al. (2017) are express as follows 


V.o+b=0 
53 
a | q) - Zo +eAs=0 G 


Ger Gerr 


2.7 Similarities and differences between gradient-damage and phase 
field methods 


In mathematics, the damage-based gradient-damage method (de Borst and Verhoosel 2016) 
and the phase-field methods for fractures are similar and these two types of methods have 
almost the same structures. Therefore, the difference between gradient-damage models and 
phase-field models is mainly in their interpretation and the interior length scale. Both damage- 
based gradient-damage and phase field methods introduce an intrinsic length scale into the 
discretization to smear the fracture over a localization band of finite width (Zhou et al. 2018c). 
However, in more detail, gradient-damage models used a spatial averaging operator and a local 
damage concept, whereas the phase field models follow the regularized energy variation due to 
fracture evolution and use a thermodynamic driving force for the smeared fracture. In addition, 
the vanishing derivative of the degradation function in the phase field ensures that the crack 
cannot broaden once the crack forms. Another difference is that crack broadening can be not 


predicted by the gradient-damage models (de Borst and Verhoosel 2016). 


3 PFMs coupled with different discretization methods 


3.1 Finite element method 


Most of the phase field methods are solved within the framework of the finite element methods 
(FEM) because the governing equations in the phase field models are commonly seen partial 
difference equations that can be solved for by using FEM. The space domain can be easily 
discretized by finite elements and suitable decoupling technology can be used to solve the fully 
coupled phase field fracture problems. Some FE implementations of the phase field models 
of fracture can be seen in Amor et al. (2009); Bourdin et al. (2012); Ehlers and Luo (2017); 
Heister et al. (2015); Hesch and Weinberg (2014); Kuhn and Müller (2008); Lee et al. (2016); 
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Liu et al. (2016); Miche et al. (2010a); Miehe and Mauthe (2016); Miehe et al. (2015b, 2010b); 
Mikelié et al. (2015a,b); Santillán et al. (2017); Sargado et al. (2017); Wheeler et al. (2014); 
Wick et al. (2016); Wu (2017); Yoshioka and Bourdin (2016); Zhou et al. (2018a,c,d). In the 
following section, more details about the implementation by using the finite element methods 
will be demonstrated. 

Although FEM has good capacity of implementing the phase field models, it is still difficult 
to solve high-order phase field equations by FEM. Therefore, some researchers tried to couple 
the phase field methods with other discretization technologies such as the isogeometric analysis 
and meshfree technology. More details about the isogeometric and meshfree methods will be 


presented in the following two subsections. 


3.2 Isogeometric method 


Borden et al. (2014) proposed a fourth-order model for the phase-field approximation of the vari- 
ational formulation for brittle fracture. Based on energy balance, thermodynamically consistent 
governing equations are derived for the fourth-order phase-field model by using the variational 
fracture approach. The higher-order model has higher regularity in solving for the exact phase 
field. To implement the higher-order model, Borden et al. (2014) employed the isogeometric 
analysis framework and the smooth spline function spaces are used. The proposed higher-order 
model improves the convergence rate of the numerical solution and complex 3D crack propaga- 
tion can be captured by the model of Borden et al. (2014). 
The modified energy functional in the model of Borden et al. (2014) is expressed as 


E(u,s) — f ((s? + qj) Vj (e) + V; (e)) ana, f 


1 1 1 
P — 8s} + ¿Vs + ze (As) dQ (54) 
Q 


4e 
By replacing Vj by the history field H, the resulting evolution equation of the higher-order 


phase field model is 
4esH 


Ge 


Schillinger et al. (2015) tried to use the isogeometric collocation methods for the discretiza- 


+s — 2@As + “A(As)=0 (55) 


tion of the second-order and fourth-order phase-field fracture models. The used isogeometric 
collocation methods are shown to speed up the phase-field fracture computations over general 
isogeometric Galerkin method because point evaluations are reduced. In addition, a hybrid 
collocation-Galerkin formulation is recommended by Schillinger et al. (2015) because it pro- 
vides consistent weakly enforcing Neumann boundary conditions and multi-patch interface con- 
straints. The hybrid collocation-Galerkin formulation can also deal with multiple boundary 
integral terms that arise from the weighted residual formulation and improve the phase field 
resolution (Schillinger et al. 2015). 
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Ambati and De Lorenzis (2016) applied a phase field model to investigate fracture in shells. 
The shell is modeled based on solid-shell kinematics with small rotations and displacements. 
The displacement and phase fields are discretized by using quadratic Non-Uniform Rational B- 
Spline basis functions. Goswami et al. (2020b) proposed a fourth-order phase field method with 
cubic-stress degradation function. Hybrid-staggered approach is used to solve the model, and 


PHT-splines are used to implement the adaptive refinement. 


3.3 Meshfree method 


Amiri et al. (2016) applied a fourth order phase-field model for fracture based on local maximum 
entropy (LME) approximants. Fourth order phase-field equations can be directly solved without 
splitting the fourth order differential equation into two second order differential equations due 
to the higher order continuity of the meshfree LME approximants. 

The higher order crack surface density function used in Amiri et al. (2016) is identical to 
that in Borden et al. (2014). The meshfree discretization form of the displacement u(a) and 


phase field s(w) is shown as follows, 


u(x) = XN paua 
s(x) = Y pasa 


where pa are LME basis functions, and ua and s, are nodal values of the displacement and 


(56) 


phase fields. The same discretization is used on virtual displacements and virtual phase field 
parameters. In addition, the LME basis functions are non-negative and must satisfy the C? and 
C! consistency: 


p(x) = 0 (57) 
Y pal) = 1 (58) 
NOD =g (59) 


where the vector £a denotes the positions of the nodes associated with each basis function. Amiri 


et al. (2016) used the local maximum entropy basis functions as follows, 


Pale) = Zu EA ea + A a — x*)| (60) 
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where 


N 
Z(a, A(z)) = Y exp [—Ay|a — z]? + A(x — a) (61) 
b=1 
is a function associated with a set of nodes X = s). y and A* is defined by 
A = in logZ(a,A 62 
arg min logZ (æ, A) (62) 


with d the dimension of the crack problem. The discrete stiffness and force matrices are then 
obtained by using the meshfree discretization to solve the crack problems. Note that the crack 
can propagate, branch and merge but cannot be recovered because a strong constraint s; < Si—1 


is imposed with s; and s;_; are the phase fields at step ? and 7 — 1. 


3.4 Physics informed neural network 


Goswami et al. (2020a) proposed the first physic informed neural network for phase field modeling 
of fracture. It seems that the PINN model treats the fracture problems one level higher on the 
model hierarchy than FEM. The proposed machine learning method is easy to implement and 
only a few lines of code are required. The approach can have a good computational savings 
compared with FEM after the network being trained. 


3.5 FFT solver for phase field modeling 


Chen et al. (2019) introduced an FFT algorithm to solve the phase-field model of brittle fracture. 
By using a staggered update scheme, the FFT algorithm solves the fracture and mechanical 
problems separately. The proposed method inherits the advantages of classical FFT methods in 
terms of simplicity of mesh generation and parallel implementation. In addition, Ma and Sun 
(2020) proposed an FFT-based solver for higher-order and multi-phase-field fracture, while the 
model was applied to strongly anisotropic brittle materials. 


4 Finite element implementation of phase field methods 


In this section, more details about the FE implementation of the phase field models are presented. 
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4.1 FE approximation 


To show better of the FE approximation, we start from the governing equations in Zhou et al. 
(2018a) for example: 


Ocij 

On; 

_ po 3 — k)H 
? Ox? G. 


+ b; = pú; 
(63) 


ES k)H 


| i Ó 
Note that these governing equations are used for dynamic problems. For quasi-static prob- 
lems, the initial term vanishes. The weak form of the governing equations are subsequently given 
by 
] (Cri su o ii ane fo duda + f:óudS = 0 (64) 
Q Q 


O» 
and 


f —2(1 — kK)H(1 — ¢)6¢dQ + f 8: (wwe . Võġ + x) dQ =0 (65) 
Q Q lo 


The standard vector-matrix notation is used and the discretization of the displacement and 


phase field can be expressed as 


u=N,d, ó-N;ó (66) 


where u; and $; are the nodal values of the displacement and phase field. d and $ are the 


vectors consisting of node values u; and ¢;. N, and Ng are shape function matrices: 


N 0 0... N 0 0 
N.=| 0 N 0 0 N, ol, Ne=| M No ... Na (67) 
0 0 M 0 0 M 


where n is the node number in one element and N; is the corresponding shape function. The 


same discretization is assumed on the test functions: 
du = N,ód, ôb = N09 (68) 


where ôd and dd are the vectors consisting of node values of the test functions. 
These gradients are used: 


e= B,d, Vó- Bb, 6¢=B,6d, Vó — Bid (69) 
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where B,, and By are the derivatives of the shape functions defined by 


Me 0 © du Me 0 0 
0 N Do... 0 N, 0 
0 o N. 0 o N, Mus Nas Nag 
B, = i oes il , B¿= Niy Noy naa Nay (70) 
Nus Nis O wee Mo Ma i 
Ni Noz Nag 
© Me Dip a 0 Mie Dag 
Ne 0 Nue Naz O Nos 


The equations of weak form (64) and (65) are subsequently written as 


—(dd)* | f pNiN,dQd + i BID, Bul + (9d)* | f Nbd + | NI fas =0 (71) 
Q Q Q 


Qh 


; e ; ; 
-00 f (BIGAB, + N¿ E + 2(1= pa No) dO + cy [ 2(1—k)HN¿d0 = 0 
Q 0 Q 
(72) 
Because of admissible arbitrary test functions, Eqs. (71) and (72) results in the discretized 


weak form equations: 


E f pNINdOd- | BID,.B,dQd- f Nbd + | Nifds =0 (73) 
Q Q Q Qn, 
——$—_ amam 
Fire=Md Fi" =Kud Feet 
Go A 
-f (BIGAB, + Nj E + 2(1 — em No) d0p + f 2(1— k) HN; dQ = 0 (74) 
9 0 Q 
ees lf SS 
Fit=K yo Fy 


where F/"*, F/"', and Fe are the inertial, internal, and external forces for the displacement 


field and F, D" and F Pus are the internal and external force terms of the phase field (Zhou et al. 


2018a). Additionally, the mass and stiffness matrices read 
M = f pNINdO 
Q 
K, = | Bi D.B,AaQ (75) 
Q 
G 


K, = i (52G.l0B, TNI E +2(1 — pa No) dQ 
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4.2 Solution schemes 


There are two basic approaches to solve the coupling discretized equations (73) and (74): mono- 
lithic and staggered schemes. In a monolithic scheme, the displacement and phase field are 
solved simultaneously. Whereas, the phase field and displacement are solved independently for 
a staggered scheme. At a given time step, the displacement is solved first at a fixed phase field. 


By using the updated displacement, the phase-field equation is subsequently solved. 


4.2.1 Monolithic scheme 


For a monolithic scheme, the implicit Hughes-Hibert-Taylor (HHT) time integration can be 
adopted in the dynamic calculations (see Liu et al. (2016)). According to our experience, 
other integration methods such as the implicit BDF (backward differentiation formulas) and 
the generalized-a@ methods can also be successfully used in phase field modeling. As an example, 
the implicit HHT method is presented herein. Letting ? be the Newton iteration, the residual 


force vectors now are re-written as (Liu et al. 2016): 
(=) Ga oa RO ate (76) 


(ROA = UR — (1 + a) (RM) + aC? (77) 


where a € [3, 0] is the HHT integration operator and At = t,, 41 — tn is the time step size. If the 
material damping effect is neglected, the tangent matrix in one element is given by (Liu et al. 


2016): 
dú 


dà — 1 _ (17 ay? 
du BAR? P5 (79) 


4.2.2 Staggered scheme 


In a staggered scheme, greater flexibility exists in solving the displacement and phase field. 
This means both implicit and explicit time schemes can be used. In an implicit approach, the 
procedure of solving displacement field is similar to those in a monolithic scheme. Whereas for an 
explicit method, e.g. the explicit central-difference time-integration method, some modification 
is needed (Liu et al. 2016): 


(u); t^ = (u); + (u); 9^! AL (80) 
(posat (ay O84" y (fA (81) 

and 
(à); = (Mi) (EE): — (Fih) (82) 
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Seles et al. (2019) discussed the details of a residual control staggered solution scheme for the 
phase-field modeling of brittle fracture. The use of a stopping criterion within the broadly used 
staggered algorithm was discussed. The methods of stopping criterion of the iterative scheme are 
based on the residual norm and implemented in commercial software ABAQUS. The proposed 
method can avoid fine loading incrementation to produce an accurate solution in a common 


single iteration staggered algorithm. 


4.3 Implementation codes 


To reduce the implementation effort on the phase field methods and facilitate new learners, some 
phase field models are developed and implemented with commonly used commercial software. 


Commonly seen are the Abaqus, FEniCS and COMSOL implementations. 


4.3.1 Abaqus implementation 


Liu et al. (2016); Msekh et al. (2015) implement the isotropic and anisotropic types of phase field 
methods in Abaqus. For the monolithic scheme, a UEL subroutine is used. The time integration 
and element discretization are implemented in the Abaqus standard and the UEL subroutine 
calls for each iteration in a given increment. Because Abaqus itself does not support a plot of 
the results on the used elements, external visualization software is required for better exporting 
the calculated results (Liu et al. 2016). 

For a staggered scheme, A UMAT or VUMAT subroutine is used. If an implicit time inte- 
gration scheme is used, a UMAT subroutine is employed to exchange the local history field H 
and the phase field $ on integration points. Subsequently, the stress ø and its tangent modulus 
with e must be renewed in every iteration. Whereas, if an explicit time integration scheme is 
used, a VUMAT subroutine is employed to exchange the local history field H and the phase 
field ¢ instead, and only the stress ø needs to be renewed during each increment. 

Zhang et al. (2018) implemented an iterative scheme in Abaqus for cohesive fractures. Fang 
et al. (2019) presented Abaqus implementation procedures of phase field fracture of elasto- 
plastic solids in a staggered manner. The subroutines UEL and UMAT are also used. The 
UMAT describes the constitutive behavior of elasto-plastic solids, while the UEL is designed for 
the phase field fracture. The authors solved the phase field and displacement field separately 


using the Newton-Raphson iteration method. 


4.3.2 COMSOL implementation 


Zhou et al. (2018a,c,d) implemented the anisotropic type of phase field models in the multi- 
field simulator-COMSOL Multiphysics. The implementations in COMSOL are much easier to 


be extended to problems with more fields than other software such as Abaqus. The modules 
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established in COMSOL is shown in Fig. 2 and the implementation procedure in COMSOL is 
shown in Fig. 3. 

As shown in Fig. 2, the “Storage Module" stores the results obtained from the “Solid Mechan- 
ics Module", such as the principal strains and their directions. The positive part of the elastic 
energy is then calculated to update the local history strain field. Subsequently, the updated 
history strain is used to solve the phase-field. Because of the highly nonlinearity, the elasticity 
matrix in the “Solid Mechanics Module" must be modified by the updated phase-field solu- 
tion and varying principal strains and their corresponding directions. Note that only staggered 
schemes are used in Zhou et al. (2018a,c,d). In addition, the implicit Generalized-a method is 
used for time integration. For convergence issues, Anderson acceleration technology is used to 
increase the convergence rate in COMSOL and the open-access codes of the phase field modeling 
in COMSOL can be found in “https: //sourceforge.net /projects/phasefieldmodelingcomsol/ ". 


Phase y € Principal strain 
> pete? 9 Direction of principal strain 
LN 


Solid. 
History field H Storage 
MI UM 
a Positive part of the 
elastic energy Components of D ] 


Fig. 2 Established modules in COMSOL for the phase field modeling 
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Fig. 3 Implementation procedure in COMSOL for the phase field modeling (Zhou et al. 20182) 
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4.3.3 FEniCS implementation 


Natarajan et al. (2019) presented a FEniCS implementation scheme of the phase field method 
for quasi-static brittle fracture. The FEniCS provides a framework for the automated solutions 
of the partial differential equations. Therefore, the phase field and displacement field can be 
solved, respectively. However, in recent days, the performance of FEniCS implementation has 


not been verified in comparison with Abaqus and COMSOL implementations. 


4.4 Computational cost 


To obtain accurate results predicted from the phase field models of fracture, it is heavily time- 
costing even in 2D problems when the finite element method is used. One reason is that the 
resolution of the small length-scale requires extremely fine meshes, at least near the region with 
high phase field ¢. That is, if no adaptive mesh refinement technology is used, the phase field 
models are calculated on fine fixed meshes, which requires high computational cost. Another 
reason is that energy split technology is used to obtain realistic crack patters under tension and 
compression. The special split such as the strain decomposition (Miehe et al. 2010a) results in 
a highly non-linear constitutive model of stress to strain. Rather many iterations are needed to 
solve the highly non-linear constitutive model. In addition, to obtain a better crack resolution, 
the loading and displacement increments must be small enough. Too small increments make the 


computational cost expensive even for an simple isotropic phase field model. 


4.5 Element technologies 


The phase field methods are in general laborious approaches to fracture and thereby some 
element technologies are developed to improve the efficiency of the phase field methods within 
the framework of FEM. 


4.5.1 Adaptive mesh refinement 


To develop a robust and efficient numerical scheme for the phase field modeling, a primal-dual 
active set method and predictor-corrector mesh adaptivity technology were used in Heister et al. 
(2015). The primal-dual active set method is a semi-smooth Newton method, which maintains 
the crack irreversibility as a constraint and achieved increasing converging rate. The predictor- 
corrector mesh adaptivity technology is also developed by Heister et al. (2015) to reduce the 
computational cost. This adaptive mesh refinement technology has following features: 

e Keep a single fixed, small strain e during the entire computation. 


e Ensure h < e inside crack region. 


e Error is controlled by e, not the element size A. 
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e No requirement of prior knowledge about crack paths. 


e Handling fast growing cracks. 


Wick (2016) developed a posteriori error estimation and goal-oriented adaptive mesh refine- 
ment technology for phase field crack propagation. Goal functionals and dual-weighted residual 
(DWR) methods are used simultaneously. Their approach is based on a partition-of-unity (PU) 
and does not require strong residuals nor jumps over element edges. 

In addition, Badnava et al. (2018); Zhou and Zhuang (2018) used another adaptive mesh 
refinement method to couple the phase field method (so-called h-adaptive phase field method). 
The adaptive mesh refinement does not require the mesh refinement around a pre-fixed crack 
path. Instead, a predictor-corrector scheme for mesh adaptivity is used. First, the coupled 
system with coarser elements is solved and predicts an initial crack path. Subsequently, a 
refinement threshold such as ¢ > 0.85 is used to judge if the mesh requires refinement. Finally, 
all the fields are solved again based on the refined mesh (Zhou and Zhuang 2018). The predictor- 
corrector mesh refinement is continuously imposed on the elements in one time step unless the 
elements do not need refinement. Tian et al. (2019) proposed a novel hybrid adaptive finite 
element phase-field method (ha-PFM) for fractures under quasi-static and dynamic loading. 
The method refines adaptively the meshes based on a crack tip identification strategy while the 
refined meshes in the noncrack progression region are reset as coarse meshes. The proposed 
method prominently reduces the CPU time and memory usage. Jansari et al. (2019) used a 
recovery based error indicator and quadtree decomposition to establish an adaptive phase field 
method for quasi-static brittle fracture. Kristensen and Martínez-Pañeda (2020) investigated the 
potential of quasi-Newton methods in facilitating convergence of monolithic solution schemes 
for phase field modelling. A new adaptive time increment scheme is proposed to reduce the 
computational cost. The study indicates computation times can be reduced by several orders 
of magnitude. On the other hand, the number of load increments required by the staggered 
solution will be up to 3000 times higher. Noii et al. (2020) proposed an adaptive global-local 


approach for phase-field modeling of anisotropic brittle fracture. 


4.5.2 Multi-scale phase field method 


Patil et al. (2018) coupled a multi-scale finite element method (MsFEM) with the hybrid type 
of phase field model for brittle fracture problems. This multi-scale method is also coupled 
with adaptive remeshing technology and named the adaptive multi-scale phase field method 
(AMPFM). An important feature of this phase field method is that the degrees of freedom of 
coarse-mesh and fine-mesh are linked together using multi-scale basis functions during mesh 
refinement (Patil et al. 2018). The crack propagation path is automatically tracked and refined 
around the crack by using the current phase field and its increment. Another benefit of AMPFM 
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is that memory and CPU time are dramatically reduced when AMPFM is used to simulate the 
brittle fracture in heterogeneous structures with uniformly distributed small-size discontinuities. 
Gerasimov et al. (2018) adopted non-intrusive global/local approaches when the fracture was 


modeled by using the phase-field framework. 


4.5.3 New element shape function 


In another effort to reduce computational cost of the phase field, Kuhn and Miller (2011) 
developed special engineered FEM shape functions to discretize the phase field. The shape 


functions have an exponential nature and their forms are as follows, 


exp(—ó(1 + €)/4) — 1 


Na cde exp(9/2) — 1 


(83) 


exp(—(1 + £)/4) — 1 
exp(9/2) — 1 
where d is the ration between the element size and the scale parameter of the phase field method. 

Note that Eqs. (83) and (84) are only available in 1D 2-node element. For 2D, the shape 
functions for a 4-node element are obtained as the tensor products of the above Eqs. (83) 
and (84) with some changes. The results of Kuhn and Müller (2011) showed that accurate 


prediction of the surface energy can be achieved in a much lower element refinement level by 


N$(6,0) = (84) 


using the special shape functions. A drawback of the approach of Kuhn and Müller (2011) is 
that some prior information about the fracture direction is needed to construct the exponential 
shape functions. 


4.5.4 Virtual element and smooth finite element 


Aldakheel et al. (2018) proposed an efficient low order virtual element scheme for phase-field 
modeling of brittle fracture. The virtual element formulation has flexible choice of node number 
in an element that can be changed in a simulation. The potential density is written in terms of 
polynomial functions rather than the unknown shape functions for complicated element geome- 
tries. Bhowmick and Liu (2018) developed a phase field model in the framework of cell-based 
smoothed finite element method (CS-FEM). The CS-FEM is softer than the standard FEM, 
and CS-FEM is used to solve the equations that govern the continuum mechanics of solids. The 


computational cost of CS-FEM is slightly lower than the finite element counterpart. 
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4.6 Special treatments 
4.6.1 Modelling pre-existing cracks 


How to model pre-existing cracks is another important issue in phase field modelling. A straight- 
forward manner is to model initial cracks as discrete cracks using double nodes in the element 
mesh. However, the phase field can also be used to model the initial pre-existing cracks. The first 
way is to constrain the phase field by the Dirichlet boundary condition o = 1 at the pre-existing 
crack. The nodes with o = 1 can be placed on a single line or a small region. But the latter 
produces larger error on calculating the crack surface energy. Another way to model pre-existing 
crack is the approach proposed by Borden et al. (2012) where an initial history energy field is 
introduced. The history field has the following form: 


G. ( _ Ue, 2) d(a,l) < e 


Ho(x) = B de P 
0 d(x, l) >€ 


(85) 


where B = 1 x 10? and d(a, 1) is the closest distance from æ to the line l. 


4.6.2 Determining crack tip position 


Because the phase field method applies diffusive description of the sharp crack, there is no obvious 
sharp crack tip in the modelling. But for post-processing, the crack tip position is sometimes 
needed especially when the crack velocity and accumulated crack length are calculated. For 2D 
problems, the crack tip is commonly defined by using the iso-curves of the phase-field (Borden 
et al. 2012; Liu et al. 2016; Zhou et al. 2018a). For example, Borden et al. (2012) used the 
iso-curves of = 0.75 to fix the crack tip position. In addition, Nguyena and Wub adopted 
another approach. They selected all nodes that have a phase field larger than 0.8, while the 
crack tip is fixed on the node that has the maximum horizontal coordinate among those selected 
nodes. 


4.6.3 Hierarchical meshes 


Goswami et al. (2019) proposed a novel dual-mesh based adaptive phase field method for solving 
fracture problems. The implementation of the phase field model is based two sets of meshes with 
different characteristic element sizes. A coarser mesh for the displacement field and a finer mesh 
is for the phase field. To facilitate the exchange of information between the meshes, the authors 


also proposed an efficient data transfer algorithm. 
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5 Extensions and applications of the PFMs 


5.1 Ductile fracture 


Extensions of the phase field methods to ductile fractures are exclusively based on the variational 
approach to fracture. These extensions are mostly based on the coupling of phase field methods 
with models of elasto-plasticity (Miche et al. 2016). Duda et al. (2015) examined a series of 
brittle fractures in elastic-plastic solids. Other variational approaches of modelling combined 
brittle-ductile fractures can be seen in Alessi et al. (2015); Ulmer et al. (2013). The model of 
Ambati et al. (2015a) suggested a characteristic degradation function. By using this function, 
the damage is coupled to plasticity in a multiplicative format (Miehe et al. 2016). However, 
combination of the local plasticity models and the gradient-damage-type phase field modeling of 
fracture in Ambati et al. (2015a) do not meet the constraints on the plastic and damage length 
scales (Miehe et al. 2016). The model of Ambati et al. (2015a) also lacks a canonical structure 
based on variational principles. This drawback is overcome by the following work of Miche et al. 
(2015a), which couples gradient plasticity to gradient damage at finite strains. 

Miehe et al. (2016) proposed a consistent variational framework for the phase field modelling 
of ductile fractures in elastic-plastic solids at large strains. The bases of the model in Miehe 
et al. (2016) are the formulation of variational gradient plasticity in Miehe et al. (2014) and the 
original phase field model for brittle fractures in Miehe et al. (2010a,b). Ambati and De Lorenzis 
(2016) applied a phase field model to investigate fracture in shells. The solid-shell formulation 
is distinguished between elastic and plastic materials. A brittle phase-field model is used for 
elastic materials, while a ductile fracture model for elasto-plastic materials with J2 plasticity 
and isotropic hardening. 

In Ambati and De Lorenzis (2016), the total free energy functional E, for the ductile fracture 


is the sum of elastic, plastic, and fracture energy contributions: 


Bilecen hes) = | [ois p) (5) + Vee) (8) (5, Vs)] do (86) 
Q 

where V, (A) is the plastic strain energy density function with isotropic hardening variable h. €e 

and e, are the respective elastic and the plastic strain tensors with the total strain e = Ee + €). 

The similar decomposition of the energy to Eq. (17) is used 


Tije 


e 


Wt (€) = ; Ks (ir(e)? 


€ 


where K, = A + 2/3 is the bulk modulus of the material. 
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The stress degradation function g(s, p) is chosen as follows: 


g(s.p) = S" +n (88) 


where m and p are two parameters with 


E e 2[' —À 
p=, e#(t) = "E / Jee, (89) 
En crit 3 0 
and &75,,, is a threshold value, and ef? is the von Mises equivalent plastic strain. 


In addition, Borden et al. (2016) presented a phase field formulation for fractures in ductile 
materials. They introduced a cubic degradation function, which produces a more accurate 
stress-strain response prior to crack initiation. A microforce-driven governing equation is used 
to replace the general energy potential for finite deformation problems. Borden et al. (2016) 
also introduced a yield surface degradation function that accounts for plastic softening and 
non-physical elastic deformations after crack initiation. 

Shanthraj et al. (2016) applied a obstacle phase field energy model to formulate the fracture 
behavior in a finite strain elasto-viscoplastic material. The obstacle energy model can produce 
physically realistic fracture behaviors at the vicinity of the crack tip. The resulting variational 
inequality is discretized by a finite element method, and is solved by using a reduced space New- 
ton method. Shanthraj et al. (2016) also emphasized a significant decrease in the computational 
cost by using their method. Alessi et al. (2018) also proposed a phase field model for plasticity. 
Based on a minimization algorithm, the coupled elasto-damage-plasticity can be solved by using 
the proposed method. Dittmann et al. (2018) proposed a higher order phase-field model for non- 
linear ductile fracture problems. The approach can easily account for the entire range of ductile 
fracture in the framework of non-linear elastoplasticity. A novel multiplicative triple split of the 


deformation gradient and a novel critical fracture energy are involved in the proposed model. 


5.2 Cohesive fracture 


For cohesive fracture, Verhoosel and de Borst (2013) developed a phase field model for straight 
crack propagation and the numerical modeling is implemented within the framework of the 
finite element method. In addition to the displacement and phase field fields in a general phase 
field method for brittle fractures, an auxiliary field was used to represent the displacement jump 
across the crack. This third field is kept constant orthogonal to the crack. Vignollet et al. (2014) 
extended the work of Verhoosel and de Borst (2013) and an arc length method and staggered 
scheme were used. However, the model of Vignollet et al. (2014) requires a pre-defined path for 


the fracture propagation. 
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5.3 Dynamic fracture 


Numerical modelling of dynamic fractures in solids based on sharp crack discontinuities (e.g. 
XFEM) is difficult in dealing with complex crack topologies and requires special branching 
criteria such as XFEM. The phase field based dynamic modelling can overcome these drawbacks. 
Following the static phase field model of Miehe et al. (2010a,b), Hofacker and Miehe (2012, 2013) 
proposed a computational framework of phase field modelling of diffusive dynamic fractures. The 
proposed modelling method allows complex crack patterns. The dynamic approach follows the 
history energy field introduced by Miehe et al. (2010a). This auxiliary field contains a maximum 
reference energy obtained in the deformation history and drives the diffusive crack evolution. In 
addition, Hofacker and Miehe (2012, 2013) applied the energy split scheme introduced in Miehe 
et al. (2010a). Some representative examples of complex crack patterns are also presented by 
Hofacker and Miche (2012, 2013). 

Borden et al. (2012) proposed a compact phase field description of dynamic fractures. The 
presented method is similar to that proposed by Hofacker and Miehe (2012, 2013). Both mono- 
lithic and staggered time integration schemes are presented by Borden et al. (2012). For the 


dynamic description, the kinetic energy V,;, is considered: 


1 


where ù = ou and p is the mass density of the material. 


The final energy functional is modified as 
T (1— sy? Os Os 
REZI | 4e noe) py 


V, = [(1— k)s? +k] UF + V; (92) 


with 


The resulting governing equations of the dynamic problem are expressed as follows, 


V-o = pù 
Ae(1 — k)H 
Go 


, (93) 
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Nguyen and Wu (2018) extended the phase-field cohesive zone model for static fracture to dy- 
namic fracture in brittle and quasi-brittle solids. The performance of the dynamic model is tested 
on several benchmarks for both dynamic and cohesive brittle fractures. Geelen et al. (2019) also 
extended the phase field formulation to dynamic cohesive fracture. The model is character- 


ized by a regularized fracture energy that is linear in the phase field, as well as non-polynomial 
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degradation functions. The authors examined two categories of degradation functions. Ren et al. 
(2019) proposed an explicit phase field model for dynamic brittle fracture. The mechanical field 
is integrated with a Verlet-velocity scheme, and the phase field is incremented with sub-steps at 
each step. Adaptive sub-stepping is applied by using the phase field residual, and the explicit 
scheme avoids the numerical difficulty in convergence and the calculation of anisotropic stiff- 
ness tensor. In addition, the phase field modulus is used rather than conventional phase field 


viscosity. 


5.4 Finite deformation fracture 


The variational approach for brittle fracture such as the original regularized formulation of 
Bourdin et al. (2000) can be extended to the problems of finite deformation fracture. Kuhn 
(2013) modified the original phase field model to describe fractures in Neo-Hookean materials 


and the elastic energy density becomes: 


A A 
V, = (s? + y) É (J?-1) — E + u)JinJ + 5 (tr(C — 3) (94) 
where J is the determinant of the deformation gradient, and C is the right Cauchy-Green tensor. 
If the surface energy density y(s, Vs) is unaltered, the thermodynamical restriction yields the 


constitutive relation: 


OW, 


= y 
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À 
= (s? +7) aun —1)C^!-4 u(1 — e» (95) 
where S is the second Piola-Kirchhoff stress tensor, and S = F~!P with P being the first Piola- 
Kirchhoff stress tensor. Being constructed in the reference configuration, the balance equation 


of the phase field modelling reads 
V-P+f=0 (96) 


where fo is the body force vector. 


Finally, the modified evolution equation of the crack phase field is expressed as 


P l8 
mc Euer 1 = — Ge 2 
M^ 312 (J ) - (à + 2u)InJ + utr(C) s)| G (eas + P ) (97) 

Note that all the parameters used in Eq. (97) are consistent with those in the model of 
Kuhn and Müller (2008). Hesch and Weinberg (2014) established another phase-field method 
for finite deformations and general nonlinear materials. A multiplicative split of the principal 
stretches is used to account for the fracture behaviors under tension and compression. Hesch 


and Weinberg (2014) also used an energy-momentum consistent integrator and their phase field 
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model is thermodynamically consistent. 


5.5 Anisotropic fracture 


Anisotropy is inherent to crystalline materials (among others) due to the symmetry of the atomic 
lattice (Nguyen et al. 2017). Failure anisotropy seemingly conflicts with the local symmetry and 
maximum energy release rate used in commonly used phase field model such as the model of 
Miehe et al. (2010a,b); Zhou et al. (2019a). Therefore, Li et al. (2015) proposed a variational 
phase-field model for strongly anisotropic fracture. In their model, higher-order phase-field 
description is implemented in a direct Galerkin way with smooth local maximum entropy ap- 
proximants. 

Nguyen et al. (2017) proposed a phase field model that could simulate non-free anisotropic 
crack bifurcations in a robust and fast implementation framework. The key in the model of 


Nguyen et al. (2017) is an introduction of the new crack surface density function: 


1 l 
10, Vo, w) = 5-0 + gw : (Vó & Vo) (98) 


where w is a second-order structural tensor that characterizes the type of anisotropy and w is 


an invariant with respect to rotations. w can be expressed as follows, 
w = 1 + 684(1— M & M) (99) 


where M is the unit vector normal to the preferential cleavage plane, and parameter £5 > 0 is 
used to penalize the damage on planes not normal to M. In the case of isotropic material, Pa 
is naturally equal to 0. 

Bryant and Sun (2018) proposed a mixed-mode phase field fracture model in anisotropic 
rocks with consistent kinematics. The mixed-mode driving force of the phase field is obtained 
by balancing the microforce. In the method, local fracture dissipation determines the crack 
propagation and kinematics modes. Pillai et al. (2020) proposed an anisotropic cohesive phase 


field model for quasi-brittle fractures in thin fibre-reinforced composites. 


5.6 Plate and shell fractures 


In this section, we provide an overview of phase field approaches applied to dimensionally reduced 
continua according to their complexity from linear to nonlinear regimes. Generally, phase field 
models applied in solid bodies can be applied to thin structures, however, the most difficulty is 
how to make quantities in phase field modeling consistent with ones in shell kinematics which 
are often described on local frames, i.e, curvilinear coordinates. On the other hand, what makes 


it different between published works until now is related to types of shell models, of phase field 
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models and of discretization methods used, which are summarized in Table 2. Here, we adopt 
the notion that classifies types of phase field model according to Miehe et al. (2010a), in which an 
approach that splits strain energy and stress tensor into negative and positive parts is regarded 
as anisotropic phase field model. While a model whose entire elastic energy is degraded in 
fracture zones is called isotropic, one that still keeps linear in the momentum balance equation 
but performs a decomposition on the elastic strain energy is referred as a hybrid phase field model 
Ambati et al. (2015b). We also note that details of related shell models are not emphasized in 


this section, instead, we refer to other review papers of shell models, see e.g Bischoff et al. (2018). 


Table 2: Overview of phase field approaches in thin structures 


Model name Shell model Phase field model Discretization method 
Ulmer et al. (2012) Mindlin-Reissner anisotropic standard FEM 
Amiri et al. (2014) Kirchhoff-Love isotropic Local Maximum-Entropy meshfree method 
Kiendl et al. (2016) Kirchhoff-Love anisotropic NURBS-based isogeometric analysis 
Ambati and De Lorenzis (2016) Solid-shell anisotropic NURBS-based isogeometric analysis 
Areias et al. (2016c) Corotational shell hybrid standard FEM 
Reinoso et al. (2017) Solid-shell isotropic standard FEM 


The first work that applied phase field model in plates and shells is proposed by Ulmer 
et al. (2012), in which the shell model is considered as a combination of standard plates and 


membranes, which allows them to use standard Lagrangian polynomials to approximate solution 
fields as the shell model only requires CO0— continuity of basis functions. On the other hand, 
elastic energy is split into bending part and membrane part which further is divided into positive 
and negative components. Accordingly, full of bending energy and positive membrane energy 


contribute to fracture progression. The free energy functional of this approach takes the form of 


= E(u.d) = f WER (es)  &)) + Enlem) + Fd, va V, (100) 
= / 

where Em is the membrane energy, € bending energy, F(d, Vd) fracture energy, and g(d) the 
degradation function. 

In 2014, Amiri et al. (2014) applied the isotropic phase field model to Kirchhoff -Love thin 
shell, for which they use Local Maximum-Entropy meshfree approximations to fulfill the require- 
ment of C1— continuity of basis functions which arises from the appearance of second derivatives 
in the bending strain of the shell model. Kiendl et al. (2016) introduced an approach that com- 
bine the anisotropic phase field model with Kirchhoff-Love thin shell and NURBS basis functions 
are used to approximate solution fields. They also shown that the two mentioned works had 
some limitations as they may prevent crack propagation and be not realistic for fracture phe- 
nomenon at some circumstances. Amiri's approach employs the isotropic model that makes it 
only applicable to cases of tensile stress states, while Ulmer’s model leads to delay of cracking 


in cases of combination between bending and compressive membrane strains or may result in 
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evolving damage under cases of purely compressive strains . Accordingly, Eq. 100 is rewritten 
as 


Haas f WOE Ee FAY (101) 
V 
where total strain € = €,, + 0? & with €m and & as membrane strain and curvature change 
respectively. In Eq. 101, both membrane and bending strains contribute to total strain tensor 
which is split into tensile and membrane components, which prohibits the decomposition of 
elastic energy into membrane and bending terms . In contrast with the aforementioned works 
that employed the classical thin shell models, Ambati and De Lorenzis (2016) investigated the 
anisotropic phase field approach to a solid-shell element whose kinematics quantities are defined 
through the thickness instead of on the mid-surface as in the Kirchhoff-Love shell. The employed 
solid-shell formulation is rotation-free, which motivates inherenting the same nodes and degrees 
of freedom of the solid element, allowing to model fracture for elasto-plastic materials by using 
general three dimensional (3D) elasto-plastic constitutive models, i.e J2 plasticity, see e.g Ambati 
et al. (2015a) for ductile fracture model in solid. 
One thing to note is that all the phase field models in thin structures above are formulated 
in linear regime, which makes them very limited , since in reality shell structures are usually 
undergone large deformation and rotation. To overcome this, Areias et al. (2016c) developed 
a hybrid phase field approach to the so-called corotational shell at finite strain. Kiendl's work 
(Kiendl et al. 2016) which ensures the irreversibility of crack evolution in the case of elasto- 
plasticity by the local history field which is the maximum of the positive elastic energy obtained 
within each loading step, Areias's model suggests another criteria for the irreversibility of crack 
phase field as 
H =< maz(W, — W5) >+ (102) 


with W, as the plastic work and W% as the critical value of the plastic work, which enables to em- 
ploy the work of separation in the ductile fracture (Siegmund and Brocks 2000). Furthermore, to 
get consistent in the extension of infinitesimal-strain regime for the continuum phase field model 
to the case of the employed finite strain shell, a consistent updated-Lagrangian algorithm, where 
a reference configuration is not chosen as the undeformed configuration, is adopted, see Areias 
et al. (2016c) for details. Recently, Reinoso et al. (2017) introduced an approach combining 
the isotropic phase field model with a six-parameter solid-shell element which consists of three 
displacements on the mid-surface, two independent rotations and additional degree of freedom 
accounting for thickness stretch. Large deformation as well as linear and nonlinear elastic con- 
stitutive laws are considered in this work. Note that, both Areias' and Reinoso's models describe 
variations of phase field through the shell thickness that leads to a correct description of fracture 


in bending-dominated cases. This is in contrast with previous works that assume invariance of 
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the phase field variable over the shell thickness. Following Areias et al. (2016c), phase field is 


linearly interpolated between top and bottom faces of the shell as follows 
1 1 
d(0*,0? 6°) = 50 + 67) drop(O", 0°) + A A 50°) (103) 


where diop and deotoom are phase field variables at top and bottom faces respectively, while 01. 
0? and 0? are coordinates in the parametric space with (0,6?) as the in-plane directions and 6? 


as the thickness direction. 


5.7 Thermal fracture 


The phase field approach can be used for thermal fracture modelling. To study thermally induced 
fractures, the thermal strains are required in the total energy functional. At small strains, the 
total strain is the sum of elastic and thermal strains. With a thermal expansion tensor œ, the 


thermal strain is €? = 0a, and the total strain e reads 
geg*-- s (104) 


where &* is the elastic strain and 0 is the temperature variation. Only the elastic part e* of the 


strain tensor contributes to the elastic energy density Kuhn (2013): 


I v= jc +q)e® : [C : e°] = (s? lle) [Ceres] (105) 


The Cauchy stress is then expressed in terms of the total and thermal strains: 


" = 40 lee (106) 


= de 2 


The crack evolution equation in Kuhn and Múller (2008) then can be modified as 


¿=-M (st E) [C : (e - e] -Ge (zeas j =<) l (107) 


In addition to the balance equation and phase field evolution equation, the governing equation 
for the temperature field is (Kuhn 2013): 


—V - ql = pc (108) 


where c is the specific heat capacity and q^ is the heat flux. Kuhn and Müller (2009) considered 


the influence of the phase field on the heat conduction behavior and the modified Fourier's law 
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can be expressed as follows 
ql = — (B(s? +n — 1) +1) KVO (109) 


where 8 € [0,1] is a parameter that defines the influence of the crack field on the thermal 
conductivity K. If 6 = 0, the heat conduction is not affected by the phase field. If 6 = 1, the 
thermal conductivity degrades and tends to zero for s = 0. Consequently, there is no heat flux 
across a crack, i.e. cracks are isolating (Kuhn 2013). 

The thermally induced phase field patterns can be seen in Bourdin (2007); Corson et al. 
(2009). Badnava et al. (2018) made a recent contribution. They adopted the similar decompo- 
sition technology of Miche et al. (2010a) to the elastic parts of the strain: 


e, = X (e)n/9n; (110) 


The initial elastic energy density is modified as 


Di = WT + wy (111) 
with i 
Ww = JEt :C: e$ (112) 


Dittmann et al. (2020) introduced a framework to simulate porous-ductile fracture in isotropic 
thermo-elasto-plastic solids and considered large deformations. In the model, they combined a 
modified Gurson-Tvergaard-Needleman GTN-type plasticity model with a phase-field fracture 
approach. Therefore, the temperature-dependent growth of voids on micro-scale followed by 


crack initiation and propagation on macro-scale can be well modeled. 


5.8 Hydraulic fracture 


PFMs seem to be a valuable alternative for modeling hydraulic fractures (HF) because all ad- 
vantages of the phase field methods are appealing for HF. PFMs for HF have been proposed 
for instance in Bourdin et al. (2012); Ehlers and Luo (2017); Heister et al. (2015); Lee et al. 
(2016); Miehe and Mauthe (2016); Miehe et al. (2015b); Mikelié et al. (2015a,b); Santillán et al. 
(2017); Wheeler et al. (2014); Wick et al. (2016); Yoshioka and Bourdin (2016); Zhou et al. 
(2019b); Zhuang et al. (2020). Wheeler et al. (2014) rewrote the energy functional by including 
poroelastic terms and succeeded in extending the phase field model to porous media. However, 
the variation in the reservoir and fracture domains with time is treated as a moving boundary 


problem and thereby extra work is needed for implementation. Later, Mikelié et al. (2015a,b) 
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fully coupled the three fields: elasticity, phase field, and pressure. They modified the energy 
functional from their previous work, the flow in their entire system is governed by Biot equations. 
The permeability tensor was also modified to consider a higher permeability along the fracture. 
The implementation approaches of Mikelié et al. (2015a,b) were then enhanced using adaptive 
element schemes by (Heister et al. 2015; Lee et al. 2016). Wick et al. (2016); Yoshioka and Bour- 
din (2016) proposed approaches to couple the phase field model to reservoir simulators. Miehe 
et al. (Miehe and Mauthe 2016; Miehe et al. 2015b) proposed new minimization and saddle 
point principles for Darcy-Biot-type flow in fractured poroelastic media coupled with phase field 
modeling. The evolution of the phase field was driven by the effective stress in the solid skeleton 
and a stress threshold was set. Moreover, the flow in the fractures was set as Poiseuille-type by 
modeling Darcy flow with an anisotropic permeability tensor. Recently, Ehlers and Luo (2017) 
embedded a phase-field approach in the theory of porous media to model dynamic hydraulic 
fracturing. Santillan et al. (2017) proposed an immersed-fracture formulation for impermeable 
porous media. 

Lee et al. (2018b) presented a framework that couples the fluid-filled fracture propagation and 
a genetic inverse algorithm for optimizing hydraulic fracturing scenarios in porous media. Lee 
et al. (2018a) proposed an immiscible two phase flow fracture model, based on a traditional phase- 
field model for predicting fracture initiation and propagation in porous media. The multifluid 
model extends the classical flow models and nonzero capillary pressure is considered. Mikelié 
et al. (2019) studied propagation of hydraulic fractures using the fixed stress splitting method. 
The mechanical step involving displacement and phase field unknowns is studied under a given 
pressure. However, these recently developed approaches for extending PFM to HF are quite 
complicated and computationally expensive. 

Another interesting work can be seen in Zhou et al. (2018c). Biot poroelasticity theory 
is applied on the porous medium with a phase field description of the fracture behavior. An 
additional pressure-related term is added to the original energy functional presented by Miehe 
et al. (2010a). However, different from Miehe and Mauthe (2016); Miehe et al. (2015b), only 
the elastic energy drives the fracture propagation and no stress threshold is imposed in the 
formulation. In addition, the phase field is used to construct indicator functions to transit fluid 
property from the intact medium to the fully broken one. Zhou et al. (2018c) implemented the 
phase field modelling by using the aforementioned Comsol Multiphysics and a staggered scheme. 
The numerical results are also verified by analytical solutions. 


The modified energy functional in Zhou et al. (2018c) reads 


But) = f 9.) - f op-(V wan f GAS- f o-u- f Finas X18) 


where p is the fluid pressure, œ is the Biot coefficient, b is the body force, and f is the surface 
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traction on the Neumann boundaries. 


The degradation of the elastic energy is modeled by the following equation: 
V.(e) = [(1 — kK)(1— 6)? +k] Vj (&) + V; (€) (114) 


where 0 < k < 1 is a model parameter that prevents the tensile part of the elastic energy density 
from vanishing and avoids numerical singularity when the phase field ó tends to 1. The same 
decomposition as Miehe et al. (2010a) is used to obtain Wj, V5, €+, and &... 

With the first variation of the energy functional and the history reference field H(x,t) = 


n wo (e(x,7)) being used, the strong form is written as 
TE(0,t 


por 


i +b; =0 
Ons (115) 
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with the Cauchy stress tensor a?" = ø (e€) — opl. 

Recently, new PFMs are developed for hydraulic fractures. For example, Shiozawa et al. 
(2019); Zhou et al. (2020) considered the effect of stress boundary in the phase field models 
and the correct HF type under stress boundary can be well predicted. Zhou and Zhuang (2020) 


proposed a phase field model for hydraulic fracturing in transversely isotropic porous media. 


5.9 length scale insensitive phase-field model 


Wu and Nguyen (2018) extended their work in quasi-brittle failure and proposed for the first 
time a length scale insensitive phase-field model for brittle fracture. Their model involves a 
phase-field regularized cohesive zone model (CZM) with linear softening law and several optimal 
characteristic functions. As an extension of the common phase field models, the length scale 
insensitive model has both failure strength and traction-separation law that are independent of 
the internal length scale parameter. The best merit of the proposed model is that it gives length 


scale independent global responses for fracture problems. 


5.10 Rock fracture 


Choo and Sun (2018) combined a pressure-sensitive plasticity model and a phase field model to 
model fractures in geological materials. Zhou et al. (2019a) revisited the formulation of Am- 
bati et al. (2015b) and established a new driving force for the phase field evolution to consider 
compressive-shear fractures in rocks. The phase field model is established in a hybrid frame- 


work. To the authors' best knowledge, the model is the first model that can simulate well the 
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compression-induced fracture in rocks and consider the friction effect. Fei and Choo (2020) 
developed another phase field model for shear fracture in pressure-sensitive geomaterials with 
emphasis on the effect of friction on fracture evolution. Governing equations for different con- 
tact conditions are established while energy is assumed dissipated during slip. The degradation 
function and threshold energy are fully considered such that the stress responses are insensitive 
to the length scale. However, the fracture direction is not automatically determined by the evo- 
lution equation of the phase field. Instead, the direction is determined using an extra criterion. 
Wang et al. (2020) proposed a phase-field model for mixed-mode fracture based on a unified 
tensile fracture criterion. The proposed method can be also applied to rock fracture. The model 
involves an additional material parameter. However, the fracture angle is determined by the 


tensile fracture criterion. 


6 Representative numerical examples 


In this section, some representative examples are presented to show the capability of the phase 


field modeling of fracture. 


6.1 Single-phasic problem 
6.1.1 Quasi-static fracture 


1. 2D notched square plate subjected to tension Fracture patterns in a square plate 
with an initial notch subjected to static tensile loading are presented in Fig. 4 (Zhou et al. 
2018a). This benchmark test has been calculated and analyzed by Hesch and Weinberg (2014); 
Liu et al. (2016); Miehe et al. (2010a,b). The geometry and loading condition are shown in Fig. 
4a. The material parameters are: Young's modulus E = 210 GPa, Poisson’s ratio v = 0.3, and 
critical eneryg release rate Ge = 2700 J/m?. Figure 4b and c shows a horizontal cracks in the 
middle of the plate. 

The effects of different phase field models on this tension example are also tested. We show 
the results obtained from Miche et al. (2010b), Ambati et al. (2015b), Amor et al. (2009) and the 
isotropic model under length scale ly = 1.5 x 107? mm and maximum mesh size h = 7.5 x 107% 
mm. The numerical simulations indicate that all the mentioned four PFMs can obtain the same 
fracture pattern. However, the fracture initiation and propagation differ at the same vertical 
displacement as shown in Fig. 5. In addition, different PFMs show different load-displacement 
curves as shown in Fig 6. Amor et al. (2009) predicts a much higher peak load compared with 
the other methods. 
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tension 
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Me 0.5 mm > 
b) lo = 1.5 x 107? c) lo = 7.5 x 107? mm 


i 0.5 mm 
(a) Geometry and boundary 


conditions 
Fig. 4 Fracture patterns of a single-edge-notched square plate subjected to tension 
(Zhou et al. 2018a) 


) Miehe et al. (2010b) 


) isotropic PFM 


) Amor et al. (2009) 


) Ambati et al. (2015b) 


Fig. 5 Fracture patterns of a single-edge-notched square plate subjected to tension when u = 


6.2 x 107? mm 
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Fig. 6 Load-displacement curve of a single-edge-notched square plate subjected to tension for 
different PFMs 


2. 2D notched square plate subjected to shear Fracture patterns in a square plate with 
an initial notch subjected to static shear loading are presented in Fig. 7 (Zhou et al. 2018a). 
This benchmark test is also shown in Hesch and Weinberg (2014); Liu et al. (2016); Miche et al. 
(2010a,b). The geometry and loading condition are shown in Fig. 7a. Figure 7b and c shows an 
inclined crack from the right tip of the pre-existing notch. As expected, the crack has a larger 
width when lo = 1.5 x 107? mm. 


M d 


0.5 mm 0.5 mm 
(a) Geometry and boundary (b) lo = 1.5 x 107? mm (c) lg = 7.5 x 107? mm 
conditions 


Fig. 7 Fracture patterns of a single-edge-notched square plate subjected to shear 
(Zhou et al. 2018a) 


The performance of the models of Miehe et al. (2010b), Ambati et al. (2015b), Amor et al. 
(2009) and the isotropic model on the shear example is also tested. Note that the same param- 


eters in the tension example are used. Figure 8 shows the final fracture patterns obtained by 
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different phase field models. Miehe et al. (2010b) and Ambati et al. (2015b) obtain a similar 
inclined shear fracture while the isotropic model and Amor et al. (2009) achieve a horizontal 
pure mode II fracture. The load-displacement curves for shear are shown in Fig. 9. Compared 
with the isotropic model and Amor et al. (2009), the plate can sustain to a larger shear load 
when the models of Miche et al. (2010b) and Ambati et al. (2015b) are applied. 


) isotropic PFM ) Miehe et al. (2010b) 
) Ambati et al. (2015b) ) Amor et al. (2009) 


Fig. 8 Fracture patterns of a single-edge-notched square plate subjected to shear: different 
PFMs 


3. 2D notched square plate subjected to tension and shear This example shows 
fractures in a double-edge-notched plate subjected to both tension and shear loading. Figure 
10a shows the geometry and boundary conditions. The length, height and thickness of the plate 
are 200, 200 and 50 mm, respectively. Two horizontal notches are 25 mm x 5 mm. More details 
can be referred to Zhou et al. (2018a). Figure 10b to 10e shows the fracture patterns under 
different critical energy release rates, i.e. Ge = 25, 50, 75 and 100 J/m?, respectively. 


4. 2D notched semi-circular bend (NSCB) test Fracture patterns in a 2D notched semi- 
circular bend (NSCB) test are simulated by Zhou et al. (2018d) and the predicted results are 
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Fig. 9 Load-displacement curve of a single-edge-notched square plate subjected to shear for 
different PF Ms 
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Fig. 10 Fracture patterns of a notched square plate subjected to tension and shear 
(Zhou et al. 2018a) 
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presented in Fig. 11. The geometry and loading condition are shown in Fig. lla. As shown 
in Fig. 11b, the fracture initiates from the upper tip of the pre-existing notch and propagates 


along the vertical direction. This observation is in good agreement with those in experimental 


tests (Gao et al. 2015). 


32mm 


A 
19mm [ — 
1 64 mm ; 1 
(a) Geometry and boundary con- (b) Fracture pattern 
ditions 


Fig. 11 Fracture patterns of a 2D notched semi-circular bend (NSCB) test 
(Zhou et al. 2018d) 


5. Propagation of multiple echelon flaws This example shows a square plate with nine 
echelon flaws subjected to tension. The geometry and boundary conditions are shown in Fig. 
12a. As observed, these flaws have the same inclination angle of 45? and a varying length and 
spacing. The final fracture pattern is shown in Fig. 12b. Fractures from the bottom flaws are 


dominated due to the stress shielding and amplification effects from flaw interaction. 


50 mm 


(a) Geometry and boundary (b) Fracture pattern 
conditions 


Fig. 12 Fracture patterns of a square plate with nine echelon flaws subjected to tension 
(Zhou et al. 2018d) 
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6. Propagation and coalescence of twenty parallel flaws This example shows a square 
plate with twenty parallel flaws subjected to tension and the fracture patterns are simulated by 
Zhou et al. (2018d) by using the anisotropic phase field model. With the geometry and boundary 
conditions being shown in Fig. 13a, symmetric fractures are depicted in Fig. 13b. In addition, 


fractures only initiate from the upper and bottom flaws and there are no interior fractures. 


qd0o00000101 


50 mm 


TONIT ITA 


(a) Geometry and bound- (b) Fracture pattern 
ary conditions 


Fig. 13 Fracture patterns of a square plate with twenty parallel flaws subjected to tension 
(Zhou et al. 2018d) 


7. Propagation and coalescence of three parallel flaws in Brazilian discs In a Brazil- 
ian splitting test, a plane cylinder specimen is subjected to diametral compression (Zhou 2018). 
Zhou (2018) used the phase field method to investigate the fracture propagation in Brazilian 
discs with multiple pre-existing notches. The specimen and flaw dimensions can be found in 
Zhou (2018). Figure 14 shows the final crack patterns of the Brazilian discs with three vertically 
arranged notches under different notch spacing. Outer cracks initiate from the notch tips and 
propagate towards the disc ends. The crack propagation intersects with the vertical direction at 
a small angle. When the spacing S = 1 cm, only one inner crack occurs between two adjoining 
notches. However, for S = 3 cm, the two inner cracks from the inner tips of the notches coalesce. 

Figure 15 shows the final crack patterns of the Brazilian discs with three horizontal notches 
and the influence of different notch spacing. For S = 2 cm and S = 3 cm, the cracks are similar. 
The phase field modelling only shows two outer cracks that propagate towards the two ends of 
the specimen. However, for S = 1 cm, additional inner cracks initiate from the notch tips and 
evolves between two adjoining notches. 


8. Compressive-shear fracture Rock is a typical geological material and compressive-shear 
fractures can be formed in rocks during loading. Zhou et al. (2019a) proposed a phase field 


model for simulating compressive-shear fracture. Here, we compare the model of Zhou et al. 
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S=1cm S=2cm S=3cm 


Fig. 14 Fracture patterns of a Brazilian disc with three vertical parallel flaws 
(Zhou 2018) 


S=1cm S=2cm S=3cm 


Fig. 15 Fracture patterns of a Brazilian disc with three horizontal parallel flaws 
(Zhou 2018) 


50 


(2019a) and the anisotropic model of Miehe et al. (2010b). The fracture pattern and load- 
displacement curves are shown in 16 under the same elastic and fracture parameters. In Fig. 16, 
the anisotropic model of (Miehe et al. 2010b) does not have a drop stage in the load-displacement 
curve while the model of Zhou et al. (2019a) has an obvious drop stage. Moreover, only wing 
and anti-wing tensile cracks are simulated in the model of (Miehe et al. 2010b) while the model 


of Zhou et al. (2019a) predicts well the compressive-shear fracture in rocks. 
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Zhou etal. (2019a) 
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Fig. 16 Comparison of fracture pattern and load-displacement curve 


6.1.2 Dynamic fracture 


1. 2D dynamic shear loading of Kalthoff experiment The phase field method has been 
used in Kalthoff experiment (Borden et al. 2012; Liu et al. 2016; Zhou et al. 2018a) and its 
geometry and boundary condition are shown in Fig. 17a. The influence of the critical energy 
release rate G, on the phase field is shown in Fig. 17b to 17e according to the results of Zhou et al. 
(2018a). As observed, a smaller Ge produces more complex crack patterns. Crack branching 
occurs for Ge = 5 x 10? and 1 x 10* J/m?. Whereas, for a larger Ge, the phase field simulation 
only show a single crack. The distance of the crack tip from the upper boundary of the plate 
increases as Ge increases. Figure 18 shows the crack-tip velocity under different Ge. As observed, 
the maximum crack-tip velocity decreases with an increasing Ge but the time of crack initiation 


increases (Zhou et al. 2018a). 


2. 2D dynamic crack branching under tension This example is a pre-notched rectangular 


plate subjected to uniaxial traction and has been investigated by Borden et al. (2012); Liu et al. 
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Fig. 17 Phase field of dynamic shear loading tests at 90 us for different Ge 
(Zhou et al. 2018a) 
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Fig. 18 Crack-tip velocity of the dynamic shear loading example for different G, 
(Zhou et al. 2018a) 
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(2016); Zhou et al. (2018a). Geometry and boundary condition of the pre-notched plate are 
given in Fig. 19. The influences of the critical energy release rate Ge on the crack pattern and 
the crack-tip velocity are shown in Figs. 20 and 21. As observed, multiple crack branching 
occurs for different G.. In addition, the crack propagates at a larger angle with the horizontal 
after the first crack branching under a larger Ge. The maximum crack-tip velocity in Fig. 21 


decreases with the increase in Ge. 


jUitttt ttt tty 
EINEN 


100 mm 


Fig. 19 Geometry and boundary conditions for the case of dynamic crack branching 
(Zhou et al. 2018a) 


Ge = 0.5 J/m?, t = 56 us ) Ge = 1 J/m?, t = 66 us 
Ge = 5 J/m?, t = 113 us ) Ge = 10 J/m?, t = 142 us 


Fig. 20 Crack patterns of the 2D crack branching example for different Ge 
(Zhou et al. 2018a) 
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Fig. 21 Crack-tip velocity of the 2D crack branching example for different G, 
(Zhou et al. 2018a) 


6.2 Hydraulic fracture 
6.2.1 Propagation of a single crack by internal fluid injection 


This example shows propagation of one single pre-existing crack when fluid is injected. The 
crack is placed horizontally at the center of a square specimen of 4 m x 4 m. The initial length 
of the crack is 0.4 m and the displacements on the outer boundaries of the square are fixed with 
the fluid pressure p = 0. Q4 elements are used to discretize all the fields with element size h = 
2 x 10? m. Here, we used the same method and parameters in Zhou et al. (2018c) to simulate 
the crack pattern and pressure field. 

The fracture pattern at t = 8.8 s is shown in Fig. 22a. As expected, the single pre-existing 
crack propagates along the horizontal direction, in good agreement with the results of Mikelié 
et al. (2015a) and Mikelic et al. (2015). The pressure distribution at t — 8.8 s is presented in 
Fig. 22b. The pressure field is similar to the phase field in Fig. 22a. However, the region where 
the pressure concentrates is larger than the cracked region because of the radial penetration of 
the fluid is used in Zhou et al. (2018c). 


6.2.2 Two parallel propagating cracks subjected to internal fluid injection 


This example shows propagation of two parallel pre-existing cracks when fluid is injected. Ge- 
ometry and boundary conditions of this example are shown in Fig. 23a. The two pre-existing 
cracks are 1 m in length and have a spacing of 0.6 m. Q4 elements are used to discretize all the 
fields with the element size h = 2 x 107? m. The evolution of the phase field is presented in 
Fig. 23b. The spacing of the propagating cracks increases with increasing time. The pressure 
distribution is shown in Fig. 23c. The pressure is similar to the phase field but has a larger 
transition band. 
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Fig. 22 Crack pattern and fluid pressure for the example of one single pre-existing crack sub- 
jected to internal fluid injection at t = 8.8 s 


6.2.3 Three parallel propagating cracks in 2D 


Crack propagation of three propagating cracks in 2D is simulated by the phase field model. 
The geometry and boundary conditions are given in Fig. 24a. The three cracks have the same 
length of 1 m and a spacing of 1 m. Uniformly spaced Q4 elements with size h = 2 x 107? m are 
employed. The crack pattern is shown in Fig. 24b. No propagation is found for the middle crack 
and only the left and right cracks propagate. The pressure distribution is shown in Fig. 24c. 
The pressure in the middle crack is much larger than that in the left and right cracks because 


the middle crack barely propagates. 


6.2.4 Propagation of two parallel penny-shaped cracks in 3D 


This example verifies the capability of the phase field modeling for hydraulic fractures in 3D. 
Two parallel penny-shaped cracks are set in the domain of (-2 m, 2 m)’. The initial crack radius 
is 0.5 m. More details can be referred to Zhou et al. (2018c). Crack propagation patterns of the 
two parallel penny-shape cracks in 3D are shown in Fig. 25. The radii of the two cracks exceed 


0.5 m when t = 13.2 s and bowl-shaped cracks are observed when t = 13.3 s. 


6.3 Fractures in plates and shells 


To illustrate the capability of phase field modeling on capturing cracking in thin structures, 
three representative numerical examples are given with a focus on qualitatively crack paths 
formed from these examples in this subsection. The employed approach is based on Kiendl et al. 


(2016), where Kirchhoff-Love thin shell and brittle fracture are considered in linear regime. To 
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(a) Geometry and boundary condition (b) Crack pattern 


(c) Pressure field 


Fig. 23 The example of two parallel pre-existing cracks 
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(a) Geometry and boundary condition (b) Crack pattern 


(c) Pressure field 


Fig. 24 The example of three parallel pre-existing cracks in 2D 
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Fig. 25 Crack patterns for the example of two parallel penny-shaped cracks in 3D 


accommodate C0-continuity requirement from the shell theory, quadratic NURBS basis functions 


are used to approximate solution fields. 


6.3.1 Single edge notched tension test 


First, a square plate consisting of a horizontal notch is considered. The top edge of the specimen 
is applied vertical displacement. The geometric properties of the plate are depicted in Fig. 26 

The material parameters are chosen such that E = 1 GPa, p = 0.3, Ge = 2 N/mm and 
lo = 0.02 mm. 23,000 elements are used to discretize the analysis domain and refine a priori in 
expected crack zones. Displacement control method is applied to conduct the simulation with 
displacement increment Au = 107% mm for the first 200 loading steps and Au = 0.5 x 1079 
mm for the remaining steps until the specimen are broken completely. Fig. 27 shows the 
crack patterns at different loading stages. The resulting crack paths are coincident with those 
obtained from previous works of Ambati and De Lorenzis (2016) and Kiendl et al. (2016). The 
load-deflection curves for solid and Kirchhoff-Love shell elements in Fig. 28 are comparable to 


each other and in good agreement with those from Ambati and De Lorenzis (2016). 


6.3.2 Simply supported Plate 


Next, a square plate without initial cracks subjected to uniform pressure is simply supported on 
its four edges. The geometry and boundary conditions of the problem are depicted in Fig. 29 
and material parameters are chosen as E = 190 GPa, p = 0.29, Ge = 0.295 N/mm and ly = 0.02 
mm. The model is discretized with 28,000 elements and h/ly = 1.7. With such applied load, 
arc-length control is utilized as a solution strategy to track mechanical response of the structure 


when crack starts to propagate. The resulting crack patterns are depicted at various stages of 
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Fig. 26 Tension test. Geometry 


Fig. 27 Tension test. Crack patterns at various loading stages 


deformation as in Fig. 30. As can be seen, crack initiation arises at the center and then crack 
branching occurs, leading to the evolution of cracks towards the corners, which is coincident 
with observations from previous investigations (Ambati and De Lorenzis 2016; Kiendl et al. 
2016). Figure 31 shows the corresponding load-deflection curve, which is in agreement with one 
obtained by Ambati and De Lorenzis (2016). 


6.3.3 Notched cylinder with internal pressure 


Finally, a curved shell structure is considered. It consists of two notches in the axial direction, 
which are located on opposite sides and the shell is subjected to internal pressure. Geometry 
and boundary condition setups are given in Fig. 32. The material parameters are chosen as 
E = 70 GPa, p = 0.3, Ge = 1.5 N/mm and ly = 0.05 mm. A mesh with 32,000 elements, which 
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Fig. 28 Tension test. Load-deflection curve 
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Fig. 29 Square plate. Geometry and boundary condition setups 


is refined a priori in the regions where cracks are expected to propagate, is used to discretize the 
shell structure . Arc-length control is adopted for the simulation. Fig. 33 shows crack phase field 
at different stages of deformation. As expected, a straight crack propagating axially is observed, 
which is coincident with those obtained from previous investigations (Ambati and De Lorenzis 
2016; Kiendl et al. 2016). 
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Fig. 30 Simply supported plate. Crack patterns at various loading stages 
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Fig. 31 Square plate. Load-deflection curve 


Fig. 32 Notched cylinder with internal pressure. Geometry and boundary condition setups 
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Fig. 33 Notched plate with hole. Crack patterns at different loading stages 


7 Conclusions 


This work presents an overall picture and some recent progress of the phase field method in 
modeling fracture initiation and propagation. This work sums up the significant stages during 
the development of the phase field models of fracture in both physics and mechanics communities. 
The theory and the use of the phase field models are also shown to demonstrate their advantages 
of handling complex fracture problems. The phase field models use an additional scalar field 
(phase field) to represent the discrete cracks rather than introduce directly physical discontinuity. 
The fracture shape and propagation are only determined from the evolution equations of the 
phase-field. Compared with other fracture methods, the phase field model does not require 
additional work to track the fracture surfaces algorithmically and thereby less computational 
effort is needed for numerical implementation of fracture. Crack branching and merging in 
materials with arbitrary 2D and 3D geometries can be easily simulated by the phase field models. 

Although the phase field model can be coupled with different discretization methods such 
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as isogeometric analysis and mesh-free methods, most phase field models use the finite element 
discretization. This work also presents details about computer implementation of the phase field 
models coupled with finite element methods. Because of the great benefits from the phase field 
modeling, the phase field models for brittle fractures are extended for solving more complicated 
problems such as ductile fractures and multi-field problems. In addition, some representative 2D 
and 3D examples for quasi-static and dynamic fractures are presented to show practicability and 
capability of the phase field modeling. For future work, the phase field models can be applied in 
fracture modeling in caverns used for compressed air energy storage (CAES) (Zhou et al. 2015b, 
2017a), which shows strong effect of cyclic stress and temperature (Xia et al. 2015; Zhou et al. 
2015a, 2018b, 2017b). 
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